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The shift in chemical equilibria due to isotope substitution is often exploited to gain insight into 
a wide variety of chemical and physical processes. It is a purely quantum mechanical effect, which 
can be computed exactly using simulations based on the path integral formalism. Here we discuss 
how these techniques can be made dramatically more efficient, and how they ultimately outperform 
quasi-harmonic approximations to treat quantum liquids not only in terms of accuracy, but also 
in terms of computational efficiency. To achieve this goal we introduce path integral quantum 
mechanics estimators based on free energy perturbation, which enable the evaluation of isotope 
effects using only a single path integral molecular dynamics trajectory of the naturally abundant 
isotope. We use as an example the calculation of the free energy change associated with H/D and 
16 0/ 18 substitutions in liquid water, and of the fractionation of those isotopes between the liquid 
and the vapor phase. In doing so, we demonstrate and discuss quantitatively the relative benefits 
of each approach, thereby providing a set of guidelines that should facilitate the choice of the most 
appropriate method in different, commonly encountered scenarios. The efficiency of the estimators 
we introduce and the analysis that we perform should in particular facilitate accurate ab initio 
calculation of isotope effects in condensed phase systems. 



I. INTRODUCTION 

Replacing an element with one of its isotopes is one 
of the most useful tools to study chemical mechanisms 
and equilibria. Since isotopes differ only by the num- 
ber of neutrons, the potential energy surface on which 
they evolve, which arises from the electron-proton inter- 
actions, is unchanged. If the nuclei behaved as classical 
particles, moving on the Born-Oppcnhcimcr potential en- 
ergy surface, the thermodynamic properties of systems 
containing different isotopes would be identical. How- 
ever, light nuclei exhibit nuclear quantum effects such 
as zero-point energy and tunnelling, which are more pro- 
nounced for lighter isotopes. As such, isotope effects pro- 
vide an experimental window through which the influence 
of nuclear quantum effects can be probed. 

Accurate and efficient simulations of isotope effects 
would enable the investigation of many important pro- 
cesses such as acid-base chemistry, shifts in phase bound- 
aries, geological fractionation ratios and equilibrium con- 
stants of chemical reactions [3, . In addition, the shift 
in the free energy barrier along a particular reaction co- 
ordinate can be used to obtain the change in the rate of a 
chemical reaction via a quantum transition state theory 
or, when combined with a dynamics scheme such 
as centroid molecular dynamics or ring polymer molec- 
ular dynamics, an approximation to the full dynamical 
rate constant [1, 0]. Isotope effects in these processes 
show a range of interesting behaviors. For example re- 
cent work on water and other hydrogen-bonded systems 
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have demonstrated that quantum mechanical effects can 
act to strengthen short hydrogen bonds or to weaken long 
ones [sj-TToj - This concept of "competing quantum ef- 
fects" has recently been shown to explain the inversion 
in the atmospherically important fractionation of hydro- 
gen and deuterium between the liquid and vapor phases 
of water where at ambient temperatures D is favored in 
the liquid whereas at higher temperatures D is favored 
in the vapor Accurately describing the cancellation 
between the quantum effects in water requires the inclu- 
sion of anharmonic terms in the potential and therefore 
conclusions drawn from simple harmonic approximations 
can be misleading. In addition the subtle nature of this 
cancellation means that fractionation ratios can be used 
as a sensitive test of a potential energy surface's ability 
to accurately predict quantum effects. 

Modelling isotope effects requires the calculation of the 
free energy change upon isotope substitution. For some 
systems such as solids and simple molecular gases a rea- 
sonable approximation of the isotope effect can be gained 
from either a harmonic approximation or h perturbation 
[l3 - ll4| . However for liquids, which contain many anhar- 
monic modes, or when treating more strongly quantum 
mechanical atoms such as H and D, these approximations 
can be inaccurate or inefficient [111 EH ■ 

The exact calculation of isotope exchange free ener- 
gies can be achieved by using the imaginary time path 
integral formalism of quantum mechanics [16H18I ]. This 
formalism exploits the isomorphism between a quantum 
mechanical system and a classical system in an extended 
"ring polymer" phase space. As such it can be applied 
to any system for which a classical simulation can be 
performed, albeit traditionally with a significant compu- 
tational overhead. 

A series of recent developments have dramatically de- 
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creased the computational cost of path integral simula- 
tions. These have taken the form of the introduction of 
efficient integrators and thermostatting schemes to ad- 
dress non-ergodicity (l9l - [2l| . In addition, methods have 
been developed to reduce the number of force evaluations 
or points on the imaginary time path [22U27I ] . These 
developments along with improvements in the ability to 
perform computationally efficient ab initio calculations 
make on-the-fly first principles path integral molecular 
dynamics a feasible proposition (28l - [3l| . Such a combi- 
nation is natural, since ab initio calculations attempt to 
reproduce the exact Born-Oppenheimer potential energy 
surface. Whereas nuclear quantum effects are often im- 
plicitly accounted for by empirical potential energy sur- 
faces, which arc typically parameterized to reproduce ex- 
perimental properties, they are completely absent in ab 
initio calculations unless they are treated explicitly. 

In this work we introduce a series of developments to 
allow the extraction of free energy changes upon isotope 
substitution from path integral simulations. Firstly, we 
show how the free energy change can be computed ef- 
ficiently by performing a physically motivated transfor- 
mation of the thermodynamic integration variables. Sec- 
ondly, we address the failures of harmonic approxima- 
tions when applied to liquids and other systems with an- 
harmonic degrees of freedom in terms of both accuracy 
and, more remarkably, efficiency when compared to path 
integral simulations. Thirdly, we introduce two path in- 
tegral estimators, based on a free energy perturbation 
approach, that offer desirable features such as the abil- 
ity to extract isotope free energy changes from a single 
path integral trajectory of the most abundant isotope. 
We benchmark these methods against the previously in- 
troduced thermodynamic integration approach [32j using 
H/D substitution in liquid water as a prototypical model 
of a hydrogen bonded system. In doing so, we demon- 
strate the strengths and pitfalls which each possess, and 
identify which approach may be most advantageous for 
a given application. 



II. THEORY 

Let us consider the thermodynamics of substituting an 
atom X of mass m in a system $ with a different isotope 
of the same element, having mass m! . More specifically, 
we will consider the process of taking an atom in a phase, 
molecule or compound, which we will label as m X($), 
and exchanging it with an atom of the other isotope taken 
from a reservoir of non-interacting atoms, which we will 
label m 'x(oo): 



»X($) + m X(oo) ^ m X($) + m X(oo) 



(1) 



The equilibrium is controlled by the free energy change 
corresponding to this process, AA Here we will consider 
the transformation to be performed at constant volume, 
but the extension to constant pressure is straightforward. 
One can compute AA by thermodynamic integration, by 



first writing out the quantum mechanical free energy at 
inverse temperature /3 for the system $ containing an 
isotope of mass \i 
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and then differentiating with respect to mass, which 
yields an expression containing the expectation value of 
the quantum kinetic energy of the system: 
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(3) 



The thermodynamic integration for the free atoms is 
straightforward, and one is left with an expression for 
the free energy difference that reads 



AA $ = ofl lo S — 

2p m 



dfx. (4) 



Here d is the dimensionality of the problem and (T(/j,, $)) 
is the expectation value of the kinetic energy of atom X 
with mass [i in the system $. 

One observation that is immediately clear from this 
result is that if the atom behaved classically, its aver- 
age kinetic energy would yield the classical equipartition 
value of 1/2/3 per degree of freedom and the free energy 
difference in Eq. Q would amount to zero. In order 
to compute AA it is therefore necessary to evaluate the 
expectation value of the kinetic energy in a way that 
takes into account the quantum mechanical nature of the 
atoms. 



A. Discretizing the free energy change 

One can in practice discretize the integral in Eq. Q 
and compute the integrand for a few values of the iso- 
tope mass, fi. In previous studies a linear discretization 
in the isotope mass [i has been used fill |32T - [34| . However, 
one can accelerate the convergence considerably by per- 
forming a change of variables that takes into account the 
physical nature of the problem. Making such a change 
of variables acts to smooth the integrand so that a good 
approximation can be obtained with a very small number 
of integration points. 

If the mass-dependence of (T(/z)) were known, one 
could make the integrand (T(/x)) j \x constant by perform- 
ing the transformation y = where f(y) satisfies 
the differential equation 



(T(f(y)))f(y)/f(y) = C, 



(5) 



where C is any constant. Clearly (T(//)) is not known ex- 
actly beforehand, so one cannot in practice solve Eq. (|5|). 
However, one can use physical intuition to make assump- 
tions on the form of (T(/i)) and derive a transformation 
which "flattens" the integral. A particularly useful limit 
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for condensed phase systems is the limit where the prob- 
lem is nearly harmonic and strongly quantized. In that 
case it can be shown that 



(0) 



where the k^s arc the spring constants of the vari- 
ous normal modes. With this assumption one obtains 
f(y) « ME.Vh^/iCy + B^Vh) 2 , where B is an 
integration constant. In practice, by choosing B = and 



the desired change of variable is simply 



ji = 1/y 2 . Eq. @ is therefore transformed into 
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(7) 



Figure Q] shows the integrand of Eqs. ^ and ([7]) for the 
H/D substitution in liquid water. Using the coordinate 
transformation in Eq. ([7J makes the integrand nearly 
constant, which simplifies greatly the integration. Us- 
ing Eq. (gl, one would obtain AA = -94.0 ± 0.1 mcV 
when computing six integration points, and —107.5 ± 
0.1 meV where only the two extrema are used. From 
Eq. ([7]) one would instead obtain the fully converged 
value AA = —93.4 ± 0.1 meV when using six points, 
and —93.5 ± 0.1 meV when using only two. 



B. The pitfalls of harmonic approximations 

For solids and gases nuclear quantum effects are usu- 
ally accounted for using the harmonic approximation [35j . 
Within this approximation, the quantum kinetic energy 
can be estimated by first computing the normal mode 
frequencies Ui and the corresponding normalized eigen- 
vectors and then evaluating the kinetic energy for a 
given atom k using [36U39} 
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(8) 



Here Uul indicates the component of the eigenvector uW 
corresponding to atom k in the Cartesian direction a. 
In solids, one would need to integrate over the Brillouin 
zone. For gas-phase molecules one can include the (clas- 
sical) contributions from rotations and translations by 
also summing over the zero-frequency eigenvalues of the 
dynamical matrix. Including quantized rotations adds 
considerable additional complexity but is only necessary 
at very low temperatures. 

For solids and simple gas-phase molecules this is of- 
ten a good approximation which can be computed inex- 
pensively. However, the cost of computing the Hessian 
grows rapidly with the number of atoms. Moreover, in 
systems such as liquids that contain many anharmonic 
modes it makes little sense to compute the Hessian for 
finite-temperature configurations |40(. However, one can 
compute the vibrational density of states for a given atom 



FIG. 1. Panel (a) shows the integrand of Eq. Q, as com- 
puted for a P = 64 PIMD simulation of liquid water with 
one hydrogen atom being transformed into deuterium. Panel 
(b) shows that the integrand becomes nearly constant when a 
transformed integration variable is used, as in Eq. ([TJ. Lines 
are just guides for the eye, and statistical error bars are re- 
ported for each point. The dashed line shows the estimate 
from the two end points. Note the difference in the scale: in 
panel (a) the integrand varies by a factor of 2, whereas in 
panel (b) it varies by less than 3%. 

a): 
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from the Fourier transform of the velocity- velocity corre- 
lation function, 



Dk(u) oc 



(v fc (i)-v fc (0))e- iwt d* 



(9) 



where is the velocity of atom k. An approximation 
to the quantum kinetic energy can then be computed 
by normalizing -Dfc(w) such that / °° Dk(oj)duj = 3 and 
evaluating [H[ 



1 h — 
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D k (uj)— coth^-dw. (10) 



Figure [2] contrasts the values of the integrand that enter 
Eq. (UJ for both liquid water and water vapor at 300 K, as 
computed exactly by direct substitution (see Sec. IIII Bp 
and by quasi-harmonic analysis (Eq. ([8]) for the vapor 
and Eq. (jTUJ) for the liquid). For each phase, the har- 
monic approximation is close to the corresponding exact 
path integral result. There are however two important 
observations to be made. 

Firstly, one is generally interested in computing dif- 
ferences between isotope-substitution free energies such 
as isotopic fractionation ratios or acidity shifts, which 
arc highly sensitive to subtle anharmonic effects. For in- 
stance the gcochcmically important process of isotopic 
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FIG. 2. The figure shows the integrand of the transformed 
thermodynamic integration © as a function of the mass of 
the isotope. The results from fully-converged PIMD (P = 64. 
blue lines) are compared with those from a harmonic/quasi- 
harmonic approximation (red lines). Curves are shown for 
both liquid water (full lines) and an isolated molecule in the 
gas phase (dashed lines) . The harmonic approximation for the 
gas phase is based on static calculations and therefore has no 
error bars. For all other calculations lines are just guides for 
the eye, and statistical error bars are reported for each data 
point. 
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the path integral simulation - the statistical error of the 
PIMD results is 10 times smaller than that of the quasi- 
harmonic approximation. In other terms, for a given level 
of statistical accuracy, a quasi-harmonic approximation 
such as in Eq. ([TO]) is 100 times more demanding than a 
fully converged PIMD simulation, even if we considered 
a number of beads twice as large as what is generally 
deemed converged and we did not exploit the efficient es- 
timators we introduce in Section [II Dl ring-polymer con- 
traction H3 or the PIGLET method [27|. 

In summary, while harmonic approximations may be a 
viable way to crudely estimate the role of nuclear quan- 
tum effects and the change in isotope substitution free 
energies in the gas phase and in simple, harmonic solids 
by means of Eq. (|5J) they should be treated with cau- 
tion in terms of both efficiency and accuracy for liquids 
and systems with anharmonicity. In particular whenever 
something as subtle as isotope fractionation ratios are 
required or whenever the system is strongly anharmonic, 
the path integral methods we discuss in the following sec- 
tion are not only more accurate, but paradoxically less 
computationally demanding than a quasi-harmonic ap- 
proximation. 



fractionation between liquid water and its vapor corre- 
sponds to the free energy change associated with the pro- 
cess 

H 2 0(Z) + HOD(u) ^ H 2 0(v) + HOD(Z) (11) 

where I denotes the liquid phase, and v denotes the vapor 
phase. Typically this is expressed as a fractionation ratio, 

ai . v = e -«AA(l)-AA(„» (12) 

Using the data in Fig. [2] to compute the fractionation 
ratio the quasi-harmonic approximation yields a value 
of 10 3 loga/_u = 156 ± 10 compared to the exact path 
integral result of 95±1. A harmonic approximation would 
therefore substantially overpredict the amount to which 
H and D separate in the atmosphere. 

Secondly, although one might view the quasi-harmonic 
approximation in Eq. (|10[) as a computationally inexpen- 
sive approach to obtain an estimate of the quantum ki- 
netic energy of a liquid, this is not the case. The reason 
for its inefficiency is that the kinetic energy is obtained 
from the integral of the Fourier transform of the velocity 
autocorrelation function. Hence any statistical noise in 
the velocity autocorrelation function must be very small 
to avoid the appearance of large statistical errors in the 
quantum kinetic energy. For example to obtain the ex- 
act results in Fig. [5] we performed a P = 64 bead path 
integral trajectory which in the most naive implemen- 
tation is ~ 64 times more costly than a single classical 
simulation of the same length. Despite the fact that we 
computed the quasi-harmonic estimate in Eq. (|10p from 
64 independent classical trajectories of the same length - 
which thereby required the same computational effort as 



C. Imaginary time path integrals 

To move beyond the harmonic approximation to the 
kinetic energy in Eq. ^ while including the quantum na- 
ture of the nuclei one can use the imaginary time path in- 
tegral formulation of quantum mechanics. This approach 
allows one to exactly include the quantum nature of the 
nuclei for a system of N distinguishable particles at a fi- 
nite temperature. For simplicity we shall use a notation 
for a single particle of mass m in one dimension, described 
by its position q and momentum p. The generalization 
to multiple degrees of freedom is routine and has been 
widely discussed elsewhere [IH, [42| . For a system with 
a Hamiltonian of the form H(p, q) = p 2 /2m + V(q), one 
can show that the quantum mechanical partition function 
Z = Tr [e-? H ] at inverse temperature /3 is isomorphic to 
the classical partition function of an extended system - 
the so-called ring polymer: 

Z P oc J dpdqe^ ffp ( p < q ). (13) 

Here Hp is the ring polymer Hamiltonian 

P_1 2 1 
Hp ( Pl q) = Yl T~ + V{qi} + 2^ {qi ~ qi+l)2 ' (M) 

i=0 

where cyclic boundary conditions i + P = i are implied, 
and up = P//3H. Equation (fl"4"| describes P replicas 
(beads) of the physical system, where replicas with adja- 
cent indices are connected by harmonic springs to form 
a closed loop. As P — >• oo the path integral partition 
function Zp converges to the correct quantum mechani- 
cal one. 
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The momenta in the partition function of Eq. (|13[) sim- 
ply provide a sampling tool and can be integrated out an- 
alytically. Hence all the quantum mechanical information 
is encoded in the configurational part of Zp. Therefore, 
even a property which is directly related to the momen- 
tum such as the kinetic energy must be computed by 
means of a configurational average of an appropriate es- 
timator. In what follows we will use (O) to indicate 
the configurational average of the estimator O(q) over 
the canonical ensemble for a particle of mass m, 



/dqO(q). 



-f E, V(q t )+mu> 2 p ( qi -q i+1 ) 2 /2 



fdq< 



(15) 



This phase-space average can be evaluated by generating 
a sequence of configurations consistent with the appro- 
priate probability distribution, and averaging over the 
values of O corresponding to the different configurations. 
This set of configurations can be generated for instance 
by means of a Monte Carlo procedure, or by molecular 
dynamics [13, EE EE 13 • 

By differentiating the partition function with respect 
to (3, one can derive the thermodynamic kinetic energy 
estimator, 

(T TD (m)) m = / — - -^mwf, ^ ( 9i - q l+1 f \ . 

\ «=0 / rn 

(16) 

In practice however computing the kinetic energy from 
Eq. (fT6"|) is problematic, because the estimator Tr_D(m) 
exhibits a large variance that worryingly grows linearly 
with the number of beads, P. Hence the statistical error 
of a simulation of a given length grows as the number 
of beads is increased to reach convergence. It is however 
possible to derive a centroid virial estimator [IE EH , 
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(17) 



where q = J^. qi/ P is the centroid coordinate. This esti 
mator yields the same average value as that in Eq. (fl6l 
but has a variance that does not grow with P. 



D. Isotope substitution by free energy 
perturbation 

As shown in Eq. (|4|) in order to evaluate the free en- 
ergy change induced by isotope substitution, one has to 
calculate the change in the quantum kinetic energy upon 
the change in mass of one of the atoms in the system. 
The most direct approach is to perform a series of inde- 
pendent simulations with different values of the isotope 
mass and compute (Tcv) f° r each. Since in most cases 
one of the isotopes has very low natural abundance, this 
typically involves performing a simulation with a single 
particle of mass fi inside a large supercell containing sev- 
eral tens or hundreds of atoms of the naturally abun- 
dant isotope. This is far from ideal, because one has 



to perform multiple simulations, one for each value of 
/i, in order to evaluate the thermodynamic integration 
in Eq. and because statistics are accumulated for 
just one atom out of the many present in the simula- 
tion. Hence a highly appealing solution would be to com- 
pute (Tcv) ^ from a single simulation of only the most 
abundant isotope by a free energy perturbation (FEP) 
method. This approach involves evaluation of estimators 
that simultaneously compute the kinetic energy and cor- 
rect the phase-space distribution to reproduce the statis- 
tics of the isotope-substituted system. 

These estimators can be obtained by writing the 
full expressions for the expectation values (Tcv) „ and 
(Tcv) m as phase-space averages. By comparing the two 
expressions, one sees that it is possible to write a "ther- 
modynamic" FEP estimator for the kinetic energy at 
mass fj, as the ratio between two phase-space averages 
computed at mass m: 



cv) 



,td (Tcv(q)exp [-h TB (^/m;q)]) r 



(exp[-/i T DW m ;q)]) r 



(18) 



where we introduce the temperature-scaled difference 
Hamiltonian 



^TD 



(a;q) 



i 



2P 



.2 P-l 

-E 



(qi -qi+iY 



(19) 



These equations provide a recipe to compute the quan- 
tum kinetic energy for an arbitrary mass ratio a = fi/m 
out of simulations performed for a single isotope mass. 
However, statistical re-weighting comes with its own set 
of problems which must be carefully scrutinized. The dif- 
ference Hamiltonian Eq. (fT9|) appears in the exponent of 
Eq. (fT5)l . If it exhibits large fluctuations the weights of 
different configurations will vary wildly and the average 
performed over many configurations will be dominated 
by a few outliers, resulting in very poor sampling effi- 
ciency. More quantitatively, recent work has shown that 
whenever one computes re- weighted averages of the form 
(ae _/i ) / (e - ' 1 ), the squared error in the mean for an av- 
erage built out of M un-correlated samples will be of 
the order of a 2 (a)e a2 ^ /M Hence the error grows 
exponentially with the variance of the difference Hamilto- 
nian, a 2 (h). In addition re- weighted averages are affected 
by a slowly-decreasing systematic error. For re-weighted 
sampling to be efficient it is therefore necessary for the 
fluctuations of the temperature-scaled difference Hamil- 
tonian to be small and certainly not much larger than 



With this in mind Eq. (|19l) is worrisome in that it is 
closely reminiscent of the thermodynamic kinetic energy 
estimator in Eq. (|16|) which is known to exhibit a vari- 
ance which grows with the number of beads, P. Due to 
this similarity we will refer to the estimator in Eq. (|T6]) 
as the TD-FEP estimator. One may therefore wonder 
whether it is possible to derive an alternative expression 
which, just as the centroid virial kinetic energy estimator 
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in Eq. dTTJ) , yields the same average with a variance which 
does not grow with P. To this aim, we use a coordinate- 
scaling approach in the form of Yamamoto (4?| ■ In our 
case, the appropriate scaled coordinates are 



Qi( a ) =1+ -7=(& -?) 

J a 



(20) 



This scaling corresponds to a contraction/dilation of the 
ring polymer. Starting from a configuration with a ra- 
dius of gyration characteristic of the abundant isotope, it 
produces one that is more compatible with what would 
be expected for an isotope of mass [i. By performing 
such a change of variables in the configurational integral 
for (Tcv) m j and then using manipulations analogous to 
those in Ref. 143, we obtain 



(Tt 



cv) 



vSC 



{T C y (q'(/i/m))exp[-fe sc (/i/m;q)]) T 
(exp[-ft gc (/x/m;q)])_ 



(21) 

Again, the average is performed in the ensemble that 
contains only the majority isotope, but now the centroid 
virial estimator is evaluated at the scaled coordinates 
fEq . (f20|) ) and the temperature-scaled difference Hamil- 
tonian is 



p-i 



h sc (a;<D = Lj2v(q'(a))-V( qi ) 



(22) 



Unlike the TD-FEP estimator in Eq. (|T^)) . it can be seen 
that the variance of hsc does not grow with P. Due to 
the scaled coordinate transformation we will refer to this 
estimator as the SC-FEP estimator. However, unlike the 
TD-FEP estimator which can be computed for an ar- 
bitrary mass ratio in a completely inexpensive way, the 
SC-FEP estimator requires an additional P energy and 
force evaluations for each mass ratio and tagged particle. 
Such force evaluations can be made inexpensive if the 
potential consists only in a sum of pairwise or otherwise 
local contributions, because in that case only the inter- 
actions involving the tagged particle need to be updated. 
However, in other cases using the SC-FEP estimator may 
entail a significant computational overhead due to the ex- 
tra potential energy calculations. 

Because of the subtle issues connected with re- 
weighting methods, in order to assess the relative merits 
of these two FEP estimators it is important to perform 
a careful quantitative analysis. In Appcndix[X]we exam- 
ine the case of a harmonic potential, where such analysis 
can be performed analytically. Let us consider here the 
resulting expressions for the fluctuations of the difference 
Hamiltonian, under the further simplifying assumptions 
that the large- P limit is being attained and that we are 
dealing with a strongly quantized normal mode having 
frequency w max (i-c. /3^w max 3> 1). For the TD-FEP 
estimator we find that <t 2 (/itd) is proportional to P, as 
one might have expected: 



c 2 (/itd) 



(a-lY 



(23) 



In contrast, the difference Hamiltonian for the SC-FEP 
estimator yields a variance which does not grow with P, 



-1 



scj 



\ 4 



(24) 



As discussed earlier, re- weighting approaches are only re- 
liable when <J 2 {h) < 1 (46|. This imposes a limit on the 
range of values of the mass ratio, a, for which the TD- 
FEP estimator can be used. The number of beads re- 
quired to converge a path integral simulation depends on 
the maximum frequency in the system and a reasonable 
convergence criterion is P = 2/3/iu; max [22| ■ Inserting this 
condition into ([23)) yields a constraint on the values of the 
mass ratio for which the TD-FEP estimator is reliable, 



1 - 



< (Xtd 1 



(25) 



Proceeding similiarly for the SC-FEP estimator yields a 
range of acceptable mass ratios of, 



1 



1 + v/8/(/3fibw - 4) 



< oisc < 



1 



1 - ^/(/Sfiow - 4) 
(26) 

These analytic estimates allow us to predict how the two 
approaches might compare in a practical case. For ex- 
ample for liquid water at room temperature the highest 
frequency is the OH stretch which extends to frequencies 
around 3500 cm -1 . This gives j3hui max ~ 16. Inserting 
this into Eqs. ([25]) and (|26| yields acceptable mass ratios 
of 



and 



0.7 < a TD < 1.3 



0.5 < a S c < 5.5 



(27) 



(28) 



for the TD-FEP and SC-FEP estimators respectively. 
Hence the TD-FEP approach in Eq. (fT8| would not even 
allow one to perform a FEP that transforms 1 H into 
deuterium (a = 2), whereas the SC-FEP approach in 
Eq. (f2Tj) would be reliable all the way to tritium (a = 3). 
In the following section we assess the accuracy of these 
harmonic predictions as to the applicability of the FEP 
estimators and contrast them with the direct substitu- 
tion method for the calculation of free energy changes in 
liquid water. 



III. RESULTS 

To assess the relative merits of the approaches intro- 
duced in the previous section we performed isotope sub- 
stitution simulations in liquid water. Despite its appar- 
ent simplicity liquid water provides a challenging test- 
bed since it contains a wide range of frequencies and ex- 
hibits a competition between two quantum effects, one of 
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which acts to strengthen its hydrogen bonding network 
and one which weakens it [|], Q . This competition is seen 
in a variety of other hydrogen bonded systems [l(J and 
arises from the anharmonicity along the hydrogen bond- 
covalent bond coordinate 0, [ll[ . As such water contains 
many of the features one may encounter in other liquids 
or hydrogen bonded systems and has an important but 
subtle anharmonic component which must be correctly 
captured in order to accurately compute the free energy 
change upon isotope substitution. 

The free energy change we will focus on is the isotopic 
fractionation of H and D between the liquid and the gas 
phase for water at room temperature defined in Eq. (jlip 
as well as the fractionation of 16 and 18 0. The results 
are quoted as fractionation factors as defined in Eq. (fT2|) 
which correspond to the ratio of isotope in the liquid 
to that in the vapor. The fractionation of isotopes is 
a frequently used technique to understand processes as 
diverse as the water cycle in the earth's climate, biolog- 
ical function and the origins of the interstellar medium 

mm. 

As can be seen from Eqs. (@| and (TTTI) the liquid-vapor 
fractionation ratio probes the change in the quantum ki- 
netic energy upon moving a water molecule from the liq- 
uid phase to the gas phase. An increase in kinetic energy 
occurs when a quantum particle is confined by the forces 
exerted on the particle by the surrounding molecules. 
Computing fractionation ratios is particularly challeng- 
ing, since their precise value depends on subtle cancella- 
tions of different quantum effects. As such, they could 
serve as a sensitive test of how well different potentials 
reproduce this delicate balance. 

We performed our simulations using the flexible simple 
point charge q-TIP4P/F water model This model 
was chosen as it has been shown to accurately predict 
the H/D fractionation between liquid water and its vapor 
[lH as well as a variety of other properties of water such 
as its density maximum and melting point when used 
in PIMD simulations. 

We will compare three different approaches for com- 
puting the free energy change upon isotope substitu- 
ion: 1.) the previously introduced direct subtitution ap- 
proach; 2.) the TD-FEP estimator in Eq. ([18]) and 3.) 
the SC-FEP estimator in Eq. (|2T]) . 



A. Generating path integral configurations 

Evaluating the expectation values of the kinetic energy 
in Eq. (fTTJ) or the re- weighted kinetic energies in Eqs. (jTSJ) 
and (|2ip requires the generation of path integral trajec- 
tories. We will compare two approaches: an efficiently 



thermostatted implementation of PIMD [21| and the re- 
cently introduced PIGLET method (2f|. 

The difference in these approaches is in their conver- 
gence with respect to the number of beads employed. In 
conventional PIMD [l8| in order to converge a simula- 
tion for a system whose fastest vibrational mode has a 



FIG. 3. a) Expectation values for the centroid-virial kinetic 
energy estimator for hydrogen in liquid water as a function 
of the number of beads. Results are shown for both PIMD 
and PIGLET simulations, and for the tagged particle having 
mass set to that of hydrogen (juh) or deuterium (mo)- b) 
Convergence of the isotope fractionation ratio between 
the liquid and the gas phase of water for the H/D substi- 
tution. Results are reported as a function of the number of 
beads, using conventional path integral MD (blue circles) and 
PIGLET (red lozenges). Lines are just guides for the eye, and 
statistical error bars are smaller than the size of the points. 
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frequency w max , the number of beads, P, must be two or 
more times the parameter /3fiu; max . For water at room 
temperature typically P = 32 yields results with an er- 
ror of a few percent. In most simulations the computa- 
tional effort is dominated by the calculation of the forces, 
and so the total cost of a conventional PIMD simulation 
scales linearly with beads - hence making the PIMD sim- 
ulation ~ 32 times more expensive than a classical one. 
The PIMD simulations were performed using the PILE- 
G thermostatting scheme [21| that has previously been 
demonstrated to be an efficient approach to simulate liq- 
uids. Since the global thermostat minimally hinders dif- 
fusion we used a short thermostat correlation time of 

Tthermo — 25 fs. 

Recently introduced colored-noise Langcvin dynamics 
approaches can reduce the cost of path integral simu- 
lations by enforcing certain pre-determined bead-bead 
correlations to be exact in the harmonic limit |2q |. This 
dramatically reduces the number of beads needed to con- 
verge many properties, even in anharmonic systems. Of 
these colored-noise approaches the PIGLET method [27j 
is particularly relevant to the present work since it was 
designed to ensure fast convergence of the quantum ki- 
netic energy, (Tcv)u which is all that is required to com- 
pute the free energy change by explicit isotope substi- 
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tution. The fast convergence of the kinetic energy using 
PIGLET is demonstrated in Figure^. Using PIMD with 
P = 32 yields an estimate of the kinetic energy with an 
error of about 2%. PIGLET yields results with similar 
accuracy to this using only 6 beads. Figure [3Jd shows 
the fractionation ratio evaluated from Eq. (fT2j) . This 
property is the small difference of two large quantities 
and is therefore very sensitive to errors in the descrip- 
tion of anharmonicities, which is apparent in the very 
poor performance of the quasi-harmonic approximation 
(see Sec. Ill Bj ). Nevertheless, PIGLET performs remark- 
ably well: from 4 to 12 beads, the values of lna;_„ differ 
by less than 10% from the asymptotic value, an error 
which is much lower than the one stemming from the use 
of an empirical potential [ill ] . 

However, even in the harmonic limit, the PIGLET ap- 
proach does not enforce all the bead-bead correlations 
needed for the expectation values of the TD-FEP and 
SC-FEP estimators in Eq. {TH} and Eq. flU). Therefore, 
there is no guarantee that PIGLET will converge much 
faster than conventional PIMD when used together with 
one of the FEP estimators. On the other hand, the ear- 
lier PI+GLE colored noise approach [26| was shown to 
enhance convergence of the kinetic energy even though 
it was not developed to enforce all the required correla- 
tions. This suggests that PIGLET, even in its current 
form, might be applied to calculate the FEP estimators 
as we investigate in Sees. IIII CI and llll Dl 



B. Direct computation of (Tcv) M 

The direct technique to compute isotope effects in- 
volves evaluating the kinetic energy (Tcv) „ in multiple 
simulations where the mass of one of the atoms is changed 
between those of the two isotopes in order to perform the 
thermodynamic integration in Eq. Q . As demonstrated 
in Sec. Ill Al with a physically motivated change of inte- 
gration variable it is possible to accurately evaluate the 
integral in Eq. (|7|) by computing just the extremal points. 
Despite the need to replicate the system this can be much 
cheaper than simply using a quasi-harmonic approxima- 
tion from a classical trajectory. 

One reason for this is the efficiency in de-correlating 
the kinetic energy in a path integral simulation. The 
statistical uncertainty of a property computed out of 
a simulation of length i s ; m decreases in proportion to 
V i~/t S i m where t is the correlation time of that property. 
Hence, for a given length of the simulation, observables 
with short correlation times will exhibit smaller statis- 
tical error. In this regard the PILE-G thermostat we 
used for PIMD calculations and PIGLET are both very 
efficient. Our simulations show the autocorrelation time 
of the centroid-virial kinetic energy estimator to depend 
weakly on P and /x, varying between 1.5 and 2 fs for all 
the simulations we performed. We used a time-step of 
0.25 fs and hence 6 — 8 force evaluations are enough to 
obtain a configuration where the kinetic energy estimator 



FIG. 4. Average kinetic energy (blue circles) computed from 
simulations in which no hydrogen atoms have been substi- 
tuted with deuterium, in a simulation box containing 64 water 
molecules. At most one atom per molecule was substituted. 
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is uncorrelated. In addition, because of the small number 
of beads required for convergence, PIGLET simulations 
could have been performed with a 0.5 fs time-step, which 
enhances the efficiency even further. This explains why 
despite the additional cost of path integral approaches 
in terms of force evaluations they are still a cheaper ap- 
proach than using the quasi-harmonic approach in Sec. 
Ill Bl In our case we used an inexpensive empirical poten- 
tial, and therefore we ran simulations that were as long 
as 500 ps, resulting in minuscule statistical errors. If one 
was prepared to accept a ±5 uncertainty on the H/D frac- 
tionation ratio, liquid water simulations of the order of 
10 ps starting from an equilibrated configuration would 
suffice, which makes performing ab initio fractionation 
calculations for liquids highly plausible. 

Due to their naturally low abundance, fractionation ra- 
tios are typically required in the low-dilution limit with 
only trace concentrations of the minority isotopes. In 
order to simulate this limit exactly the mass has to be 
changed for a single atom in a system that contains hun- 
dreds of instances of the naturally-abundant isotope, and 
the kinetic energy can be computed for that particle only. 
However, in a homogeneous system such as water it is in- 
teresting to explore how sensitive the kinetic energy, and 
hence the computed fractionation ratio, are to multiple 
isotope substitutions in the system. The advantage of 
making these multiple isotope substitutions is that one 
can then gather statistics from several atoms at the same 
time. Momentum-dependent observables tend to be very 
local in nature, and one can hope that the error intro- 
duced by having an unnaturally high concentration of 
isotope substituted species will be small. For instance, 
when performing open-path PIMD in liquid water, one 
can simultaneously open the path for one hydrogen atom 
per water molecule without introducing significant sys- 
tematic errors [Hi], [Hj]. Figure 2] shows how {Tcv) mD 
changes when multiple atoms in a simulation of light wa- 
ter are transmuted into deuterium, with the precaution 
of never substituting two atoms on the same molecule. 
Even when half the H's in the system are converted, the 
kinetic energy of the D's changes by less than 0.3%, a dif- 
ference that is not statistically significant even with these 
relatively long calculations. When performing more de- 
manding ab initio calculations, the statistical error from 
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the short runs will almost certainly be dominant, making 
it particularly attractive to gather horizontal statistics by 
multiple substitutions. One can also observe that when 
multiple substitutions are performed the statistical un- 
certainty of (Tcv) fl — mD decreases with the square root 
of the number of substituted atoms. This implies that 
the instantaneous values of Tcv for atoms in different 
molecules are almost perfectly un-correlated - another 
indication of the local nature of the quantum kinetic 
energy. However, we stress that one should be careful 
when performing multiple substitutions, in particular for 
inhomogeneous systems, and that whenever possible one 
should benchmark this approximation against the highly- 
diluted case. 

Even exploiting our efficient transformation of the in- 
tegral which requires only simulations at the two extrema 
the need to perform multiple simulations to obtain free 
energy changes using direct substitution is undesirable. 
We now consider our FEP approaches which allow the 
free energy change to be obtained from a single trajec- 
tory of the most abundant isotope. 



C. "Thermodynamic" free energy perturbation 

The TD-FEP estimator in Eq. (fT5|) provides an appeal- 
ing approach to the computation of the kinetic energy of 
an isotope-substituted system by just post-processing the 
data from a simulation that contains only the naturally- 
abundant isotope. For a homogenous system such as 
water this "virtual substitution" can be performed one 
atom at a time for every atom in the system allowing 
full horizontal statistics to be collected. However, as dis- 
cussed in Sec. |H] the re-weighting procedure which un- 
derlies free-energy perturbation methods is potentially 
dangerous, as it can introduce systematic errors and a 
difficult-to-detect statistical inefficiency j46[. The pres- 
ence of unusual statistical behavior, dominated by out- 
liers, is demonstrated in Fig.[5j which shows the cumula- 
tive average of the kinetic energy predicted from Eq. (TTg|) 
for different mass ratios as the simulation progresses. 
For large values of fx/mn the cumulative average does 
not converge with regularity to a mean value, but ex- 
hibits a characteristic sequence of plateaus, interrupted 
by abrupt jumps when a new outlier is encountered. 

A more quantitative approach to monitor the statisti- 
cal quality of the TD-FEP estimates of {Tcv)^ D involves 
computing the variance of the temperature scaled differ- 
ence Hamiltonian in Eq. (|T9")) . As discussed in Sec. |H] 
a variance, o~ 2 {h), smaller than one is required for rc- 
wcighting to be of practical use. Figure [6] shows the vari- 
ance obtained from our water simulations for different 
combinations of the number of beads, P and mass ratio, 
a. As expected from our analysis in Appendix [A] the 
variance of the TD-FEP grows rapidly with both P and 
a so in practice this estimator becomes useless for mass 
ratios greater than about 1.2. Importantly the simulated 
behavior of the variance is in excellent agreement with 



FIG. 5. Cumulative average of (Tcv) „ , as computed for 
hydrogen isotopes in a single PIMD trajectory of liquid wa- 
ter, with P = 32. Lines, from top to bottom, correspond to 
different mass ratios ranging from ^/mn = 1 to [i/mii = 2. 
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FIG. 6. Contour plot of the variance of the difference Hamil- 
tonian for the thermodynamic free energy perturbation hm, 
as a function of the number of beads P and the mass ratio 
a = (i/niH- For a 2 (hrD) > 1 computing {Tcv)^ D becomes 
impractical. 




the analytic prediction derived for the harmonic limit in 
Eq. (|27p which predicts a maximum mass ratio of 1.3 
for water simply from the knowledge of the highest fre- 
quency present. Hence, the analytic estimates provide 
an excellent a priori way of deciding whether one should 
employ TD-FEP for a given system. 

A maximum mass ratio of 1.2 using conventional PIMD 
simulations implies that TD-FEP is not recommended for 
studying isotope substitution for hydrogen in room tem- 
perature water but this mass ratio is within the range 
needed to perform 16 0/ 18 substitution, which involves 
a mass ratio of only a = 18/16. Figure [7] shows the 
isotope fractionation ratio for 16 O and 18 O between the 
liquid and the gas phases of water at ambient conditions. 
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FIG. 7. Convergence of the isotope fractionation ratio ai- v 
between the liquid and the gas phase of water at ambient 
conditions for the 16 0/ 18 substitution. Upper panel depicts 
the convergence of the 18 O kinetic energy in the two phases, 
as computed from thermodynamic free energy perturbation, 
while the lower panel reports the fractionation ratio. Lines 
are just guides for the eye. 




The isotope substitution free energies were computed by 
TD-FEP from a single simulation of pure 16 water. De- 
spite being a heavy atom the number of beads required to 
converge the kinetic energy in the liquid and gas phase 
is ~ 32 since it is bonded to a light atom. The frac- 
tionation ratio converges with the number of beads more 
rapidly than the kinetic energy estimator does, because 
of a cancellation of errors between the liquid and the 
gas phases. We find I0 3 lnai- v = 12.5 ± 0.5, for the q- 
TIP4P/F model which is in reasonable agreement with 
the experimental value of 9.2 [53| . This is in error by 
an amount similar to a fall of 30 K in the temperature 
which is roughly the same as the error in the H/D li quid 
vapor fractionation when compared to experiment [ll| 
suggesting this model mildly over-predicts fractionation 
in water. 

As we discussed in Sec. IIII AI PIGLET does not enforce 
the bead correlations necessary to ensure that it will ac- 
celerate the convergence of the TD-FEP estimator even 
in the harmonic limit. However, there are two reasons 
why such a combination could be highly beneficial and 
hence is worth investigating. The first is simply the re- 
duction in the number of beads, which implies that fewer 
force evaluations are required to evolve the system result- 
ing in a speed-up of ~ 6 times in water. The second is 
that the variance of the difference Hamiltonian that ex- 
ponentially determines the rise in the error of TD-FEP 
depends on the number of beads (see Fig. [5]). Hence 
any reduction in the number of beads by using PIGLET 
would offer exponential benefits in the efficiency of the 
TD-FEP approach. 

Figure [5^ shows the value of the kinetic energy of 
(Tcv) mD estimated from a single simulation of pure 



FIG. 8. The figure demonstrates the performance of the TD- 
FEP estimator used together with PIGLET. Panel a) repre- 
sents (Tcv) mD as a function of the number of beads. PIMD 
results are shown for comparison, and the black line indicates 
the reference value (direct substitution, PIMD, P = 64). In 
Panel b) the isotope fractionation ratio between liquid and 
vapor is shown, for water at room temperature, as estimated 
from TD-FEP and PIGLET. The black line indicates the ref- 
erence value (direct substitution, PIMD, P = 64). 
a) 
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liquid H 2 0, using TD-FEP with PIGLET. Without 
PIGLET, at least 32 beads would be needed to converge 
the kinetic energy, which as discussed above would re- 
strict the maximum acceptable mass ratios to about 1.2. 
On the other hand, as shown in Fig.[BjD, PIGLET reaches 
convergence with about 6 beads for both the kinetic en- 
ergy and the fractionation ratio, which is just on the 
boundary of the region for which performing TD-FEP 
up to a = 2 would become impractical. The fact that 
PIGLET was not parameterized to enforce the particular 
correlation needed for TD-FEP does not appear to limit 
its utility in converging it and indeed remarkably it even 
seems to converge as rapidly as kinetic energy, which is 
rigorously enforced in the harmonic limit. Although by 
reducing the required number of beads PIGLET delays 
the failure of TD-FEP the statistical error blows up ren- 
dering the estimator unusable beyond 6-8 beads (see Fig- 
ure O. The extra window of mass ratio PIGLET opens 
is however incredibly useful since it allows H/D isotope 
substitution to be computed at room temperature and for 
ubiquitous aqueous and hydrogen-bonded systems, per- 
forming just a single trajectory and using as few as 6 
beads. 
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FIG. 9. Contour plot of the variance of the difference Hamil- 
tonian for the coordinate-scaling free energy perturbation 
hsc, as a function of the number of beads P and the mass 
ratio a = fx/mH- Throughout the range we considered, 
o 2 {h sc ) < 1. 




D. Scaled coordinates free energy perturbation 



FIG. 10. The figure demonstrates the poor convergence of 
the SC-FEP estimator when used together with PIGLET, a) 
{Tcv) m as a function of the number of beads. PIMD re- 
sults are shown for comparison, and the black line indicates 
the reference value (direct substitution, PIMD, P = 64). b) 
Isotope fractionation ratio between liquid and vapor for water 
at room temperature, estimated from SC-FEP and PIGLET. 
The black line indicates the reference value (direct substitu- 
tion, PIMD, P = 64). 
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In the previous section TD-FEP was shown to be lim- 
ited due to possessing a statistical error that grows with 
the number of beads. Based on our harmonic analysis 
the SC-FEP approach should allow one to obtain a much 
wider range of mass ratios a from a single trajectory. 
In addition the analytic predictions suggest that, unlike 
TD-FEP, the variance of the difference Hamiltonian for 
SC-FEP will show no dependence on the number of beads 
allowing the path integral simulation to be performed to 
high precision without the statistical inefficiency rising 
to unacceptable levels. 

To test these predictions Fig. |H] shows the variance of 
the difference Hamiltonian hsc for SC-FEP as a func- 
tion of P and fi/ma, over the same range of masses 
and numbers of beads we explored in Figure |5] for TD- 
FEP. In striking contrast with the case of the TD-FEP, 
f 2 (^sc) is nearly independent of P, and even up to a 
mass ratio of a = 2 is well below the threshold value 
of one. The SC-FEP estimator is therefore suitable for 
computing H/D substitution using conventional PIMD 
using cither P = 32 or 64 beads depending on the accu- 
racy required. Based on the scaling of the variance with 
temperature predicted by Eq. (|26[) one would expect this 
to be the case even at temperatures considerably lower 
than 300 K. 

Given SC-FEP 's much greater range of applicability 
than TD-FEP when used in conjunction with conven- 
tional PIMD one might wonder if SC-FEP could be re- 



duced in cost by using PIGLET which was shown to be 
highly effective in accelerating convergence of TD-FEP. 
Figure [TU] demonstrates that in this case the bead-bead 
correlations needed to converge the SC-FEP estimator 
are clearly not included in the current PIGLET parame- 
terization leading to very poor convergence. When using 
PIGLET with 6 beads, which is essentially converged for 
TD-FEP or direct substitution, the SC-FEP estimator 
gives a kinetic energy estimate for deuterium with an 
error of 100% and a fractionation ratio which does not 
even possess the correct sign. It is reassuring to note 
that, by further increasing the number of beads, conver- 
gence will be eventually reached. This underlines the 
robustness of the PIGLET approach: when a sufficiently 
large number of beads are used, the non-equilibrium na- 
ture of the colored noise naturally vanishes and these 
methods reduce to conventional PIMD. It is also worth 
mentioning that in principle one could exploit the philos- 
ophy of the PIGLET method to enforce in the harmonic 
limit the bead-bead correlations that are necessary to 
rigorously converge both TD-FEP and SC-FEP. Such a 
re-parameterization of the colored noise is a non trivial 
task and hence is left for future work. 
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E. Recommendations 

From the discussion above it is clear that each ap- 
proach for obtaining the free energy change upon isotope 
exchange has benefits and deficits. As such there can be 
no single recommendation but rather, based on our sim- 
ulations and analytic results, one can suggest the most 
appropriate approach to efficiently tackle typical simula- 
tion scenarios. 

The direct substitution approach provides a baseline 
since it has been extensively used in previous studies 
[Til HH-HH . In the present paper we demonstrate a num- 
ber of stratgies that greatly improve its efficiency. For 
a start, by performing a physically motivated smoothing 
of the integral fSection lll Ap one can usually perform the 
thermodynamic integration with only two simulations for 
the extremal values of the isotope mass. In addition the 
PIGLET approach can be used to reduce the number of 
force evaluations required for each simulation. Finally, 
our results suggest that in a homogeneous system per- 
forming multiple substitutions is a reliable way to obtain 
horizontal statistics - although this should of course be 
treated with caution and tested carefully before it is ap- 
plied to other problems. With these improvements, direct 
substitution becomes a very competitive choice even in 
comparison to approximate methods. 

The free energy perturbation approaches we introduce, 
however, are even more appealing, in that they allow one 
to obtain the isotope substitution free energy from a sin- 
gle simulation containing only the majority isotope, and 
allow one to gather horizontal statistics in an exact way. 
However, one must be careful of the statistical efficiency 
of re-weighting whenever large fluctuations occur in the 
difference Hamiltonian. In the preceding sections we have 
presented both analytical and simulation evidence to de- 
termine when FEP approaches are effective. This allows 
us to provide some simple practical guidelines: 

• For heavy isotope substitution at temperatures 
above 100K, one should perform a PIGLET sim- 
ulation containing only the most abundant isotope 
and use the TD-FEP estimator to extract the free 
energy change. 

• For H/D substitution in "biological" conditions 
(room temperature, aqueous systems) and at 
higher temperatures, one should also use TD- 
FEP and PIGLET. However, this more challenging 
regime sits on the edge of the TD-FEP estimator's 
applicability when combined with PIGLET. Hence 
TD-FEP should not be used with PIMD in this 
regime since the additional beads required for con- 
vergence cause the statistical efficiency to become 
unacceptable. 

• For isotope substitution below room temperature 
for H/D or at room temperature when making a 
substitution with a larger mass ratios than H/D 
(e.g Muonium/H or H/Tritium) a direct substitu- 



tion calculation is the best option. The efficiency of 
direct substitution can be greatly enhanced by us- 
ing our transformed integration variable, PIGLET 
and, where appropriate, using multiple substitu- 
tions. 

While the three cases above cover most practical ap- 
plications, there are circumstances in which the SC-FEP 
estimator becomes competitive. Due to the requirement 
to evaluate the potential on a scaled ring polymer its 
main advantage lies in cases where the potential can be 
evaluated efficiently upon a single particle change such 
as for empirical potentials. As we demonstrated the con- 
vergence of the SC-FEP is not improved by PIGLET and 
so PIMD must be used to simulate the trajectory. How- 
ever for empirical potentials PIMD trajectories can be 
made extremely cheap by using ring-polymer contraction 
techniques (22l. [2~I1. |4~D] . SC-FEP is also more attractive in 
problems where only a single atom in the system is sub- 
stituted such as for a particular proton in an enzyme 
active site. 

Since our analytic predictions of the fluctuations in the 
difference Hamiltonians for the FEP estimators were seen 
to be reliable when compared to liquid water simulations 
we can use them to derive simple expressions that pre- 
dict the efficiency of TD-FEP and SC-FEP relative to 
direct simulation. These efficiency measures, Etd and 
Esc gauge the cost of obtaining the free energy change 
from the TD-FEP and SC-FEP estimators respectively 
compared to obtaining it from direct substitution. If 
Etd or Esc come out much greater than one then the 
FEP approach is preferable and the FEP variant with 
the highest value should be used. If both values come 
out much less than one then direct substitution should 
be employed. Values around one suggest that either FEP 
or direct substitution are of comparable efficiency for the 
problem. For an isotope substitution with a mass ratio 
a = m'/m between the most and least abundant isotope 
at inverse temperature /3 in a system with a maximum 
vibrational frequency w max one can show that for the 
TD-FEP estimator the efficiency is 
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Here ni n t is the number of integration points one needs 
to converge the thermodynamic integration (2 or 3 when 
using Eq. (J7J), Nx is the number of equivalent atoms for 
which the substitution can be performed, and Mx the 
number of multiple substitutions one can perform in a di- 
rect calculation (Mx = 1 is the only completely safe bet). 
P is the number of beads required for convergence, which 
is of the order of 2/3fiw max with conventional PIMD and 
typically about /3fiw max /3 with PIGLET. From this equa- 
tion one can immediately see that for systems which re- 
quire large numbers of integration points, rii n t or posses a 
large number of exchangeable atoms then TD-FEP holds 
an advantage. However, at low temperature (high j3) 
or a large mass ratio or high frequencies w max the effi- 
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ciency can quickly shift in favour of direct substitution. 
One can also immediately see why PIGLET is so useful 
when used with TD-FEP since the relative efficiency is 
exponentially sensitive to the number of beads used so, 
by reducing the number of beads required by a factor of 
~ 6, PIGLET yields a e 6 ps 400 times increase in Etd- 
The efficiency of the SC-FEP estimator is, 
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Besides the quantities present in Eq. ([29]) . this expres- 
sion also contains r. the correlation time of Tcv, Ai, the 
time step of the simulation when using PIGLET, and 
K, the cost of computing the scaled-coordinates poten- 
tial relative to the cost of a simulation step (e.g. K = 1 
for ab initio, and virtually K = if the potential can 
be evaluated locally). Since SC-FEP cannot be used to- 
gether with PIGLET, the expression also contains both 
the number of beads required to converge PIGLET, Pole 
and the number of beads required to converge conven- 
tional PIMD, P. Here it can be seen that for systems 
with long correlation times r, cheap scaled-coordinate 
potential evaluations K or low numbers of exchangeable 
sites Nx, SC-FEP dominates. 



IV. CONCLUSIONS 

In this paper we have introduced a series of improve- 
ments which enhance the ease and efficiency of calculat- 
ing isotope exchange free energies. This include two new 
estimators, based on free energy perturbation techniques 
as well as a number of improvements to the existing direct 
substitution approach. We then critically discussed the 
advantages and shortcomings of these approaches predi- 
cating our reasoning on both analytic arguments and on 
extensive simulations of room-temperature water. 

The balance of pros and cons is fundamentally related 
to the statistical efficiency of sampling of the different 
techniques. This, in turns, depends on the "quantum- 
ness" of the system, as expressed by the ratio between 
the highest quantum of vibrational energy and the ther- 
mal energy, /3fiw max and on the ratio of the masses of the 
minority isotope and of the naturally abundant one. 

When this mass ratio is close to one, as is the case 
whenever one is dealing with relatively heavy elements 
such as oxygen, a "thermodynamic" FEP estimator al- 
lows one to compute inexpensively the effect of isotope 
substitution based on a single simulation which only con- 
tains the majority isotope. This both simplifies and 
makes less computationally demanding the evaluation of 
quantities of great interest for geochemistry and clima- 
tology. 

When instead the mass ratio is large, such as when 
one needs to substitute deuterium for hydrogen, or when 



the temperature is considerably below room temperature, 
the scenario is less clear-cut and one has to weight care- 
fully a number of factors. The statistical efficiency of the 
"thermodynamic" estimator degrades dramatically with 
the number of beads, so it becomes mandatory to use 
PIGLET to ensure that convergence is achieved with a 
reduced number of beads. With this caveat, H/D substi- 
tution is still manageable with TD-FEP in most systems 
of interest at room temperature. At cryogenic tempera- 
tures an alternative estimator can be used, which is based 
on a scaling of coordinates relative to the centroid of the 
path and exhibit more robust statistical properties. How- 
ever it requires multiple evaluations of the forces, and 
is therefore more computationally demanding than TD- 
FEP. For simulations using empirical potentials it is pos- 
sible to reduce this computational cost by exploiting the 
locality of the coordinate scaling procedure. For an ab 
initio simulation, however, there might be corner cases 
in which a carefully-performed direct substitution sim- 
ulation is more advantageous that either TD-FEP and 
SC-FEP. 

We provide analytical expressions that allow one to es- 
timate a priori the different factors that determine the 
sampling efficiency, so that one can make the optimal 
choice for a given system without having to perform time- 
consuming benchmarks. This analysis and the estimators 
we introduce will greatly simplify the evaluation of ob- 
servables connected with isotope substitution, that can 
be measured experimentally with great precision and pro- 
vide precious insight into a multitude of natural phenom- 
ena. Obtaining the corresponding theoretical predictions 
will be useful both to elucidate and prcdic the results of 
experiments, and to validate and benchmark empirical 
and ab initio predictions of condensed-phase systems. 
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Appendix A: Reweighting statistics for FEP 
estimators 

When developing new estimators in a path integral 
context, analysing their properties in the harmonic limit 
often allows one to obtain analytic results and useful in- 
sight into their limitations. This is particularly impor- 
tant when performing re-weighted sampling, because of 
the risk of dramatically degraded sampling statistics. In 
what follows we will consider a one-dimensional harmonic 
oscillator of frequency ui and mass m, simulated at tern- 



14 



perature Its/ (3 and using P beads. Consider for instance 
the expectation value of the difference Hamiltonian for 
the TD-FEP estimator (Q~9|). One can first transform 
into the normal- modes coordinates of the path [42j . 
and get 



TD 



(a — \)mj3 
2P 



4 (qI), 



(Al) 



One can then compute analytically the average fluctua- 
tions of the normal mode coordinates, and take the limit 
of an infinite number of beads analytically, which yields 



(h 



TD/ m = 



a-l 
2 

a-l 



P^oo 



(A2) 



„ X , X 

P coth - 

2 2 



Here we introduced the frequency of the fc-th normal 
mode of the free ring polymer, ujk = 2cop sin kir/P, and 
the shorthand x = ptvuj. Proceeding in a similar way, one 
can obtain the fluctuations of /itd, 



c 2 (/itd) 



(a — l)m(3 



2P 



2 2EN(^)J 2 = 



(«-l) 2 , ^l-2(l+^/ W 2 )' 



2 

(a-l) 2 



„ 3 , x 1 / x 2 

P - -x coth - + - —r 

4 2 2 Vsinha;/2 
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As one could have imagined, fluctuations exhibit a term 
growing linearly with P, which in the context of re- 



weighted sampling has particularly disruptive conse- 
quences. 

The procedure for the SC-FEP estimator ([21]) is anal- 
ogous. One evaluates the expectation value of the differ- 
ence Hamiltonian 



sc) 



(a- 1 - 1)t 



2P 



E- 2 <<a 



A;>0 



OL~ - 1 



fc>0 

y — — 1 (X 



V 

f-' 1 + LU 
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(A4) 
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2 2 



which does not contain an explicit P dependence in the 
asymptotic limit. This is also true of the fluctuations, 
that evaluate to 



<x 2 (/isc) 



(a" 1 - l)m/3 



2P 
2 

(a- 1 - If 



fc>0 

1 



V- 

^ (1+ujI/uj 2 ) 2 p^°° 
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X X 

— coth — , 
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x/2 



(A5) 



Eqs. (|A3[) and (| A5|) capture the essential statistical prop- 
erties of the TD-FEP and SC-FEP estimators, and can 
be profitably used to assess the range of applicability of 
the two methods by computing analytically the order of 
magnitude of the statistical errors one may expect for a 
system with known normal mode frequencies. 
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